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Computing the eigenvalues of fourth 
order Sturm-Liouville problems with 
Lie Group method 


H. Mirzaei* 


Abstract 


In this paper, we formulate the fourth order Sturm-Liouville problem 
(FSLP) as a Lie group matrix differential equation. By solving this ma- 
trix differential equation by Lie group Magnus expansion, we compute the 
eigenvalues of the FSLP. The Magnus expansion is an infinite series of mul- 
tiple integrals of Lie brackets. The approximation is, in fact, the truncation 
of Magnus expansion and a Gaussian quadrature are used to evaluate the 
integrals. Finally, some numerical examples are given. 


Keywords: Lie group method; Fourth order Sturm-Liouville problem; Mag- 
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1 Introduction 


A classical fourth order Sturm-Liouville equation, is a real fourth order linear 
differential equation of the form 


(p(x)y"(a))" — (qu(a)y"(@))’ + ga(a)y(@) = Aw(2)y(@), (1) 


where p(x), qi(Z), g2(x) and w(x) are given continuous functions on the finite 
interval [a,b]. Usually the equation (1) isaccompanied with four boundary 
conditions of the form 


Ayu(a) + Agu(a) =0, Byu(b) + Bov(b) = 0, (2) 


where u = [y, y'], v = [py”, (py”)’ +qry’], Ai, Az, Bi and Bz are real matrices 
of order 2, such that 
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A, Ai = A,AT, B,BY = BBT, 


and the matrices (A; : Az) and (B, : Bz) have rank 2, see [8,9]. The equa- 
tion (1) with boundary conditions (2) is called fourth-order Sturm-Liouville 
problem (FSLP). If A is such that the FSLP has a nontrivial solution, then 
A is called an eigenvalue and nontrivial solution for corresponding to A is 
called an eigenfunction. The set of all eigenvalues of FSLP called the spec- 
trum. We assume that the interval (a,b) is finite and coefficient functions 
r= WN, q2 are real and belong to L1(a,b). The classical results of self- 
adjoint Sturm-Liouville problem states that under this assumptions the eigen- 
values are bounded below and can be ordered as 


do < AL < Ag <-+, (3) 


where limg_.o0 Ax = 00, see [8-10]. 


FSLP arise in mathematical modeling of motion for the transversely vi- 
brating Beam. In fact, the equation (1) is the standard form of Euler- 
Bernoulli Beam equation, see [7]. Moan [15] approximates the eigenvalues 
of second order Sturm-Liouville problems using Lie group methods. Ledoux 
and et al. [17] introduced a modified Magnus method for numerical solution of 
one dimensional Schrodinger eigenvalue problem. In [1] the Homotopy anal- 
ysis method (HAM) is applied to numerically approximate the eigenvalues of 
the second and fourth order Sturm-Liouville problem. For further schemes 
to study FSLP, we refer to [2,5, 10]. 


Let uw = y, ua = y’, ug = py”, us = (py”)’ — qy’. Then we have the 
matrix representation of equation (1) and boundary conditions (2): 


U! =G(a)U 
{ AU(a) + BU(b) = 0, (4) 


where 
0 1 0 O 
G(z) = 0 0 r(x) 0 (5) 
0 q(x) O 1 
Aw(x)—go(x) 0 0 0 


The paper is organized as follows. In Section 2, we introduce Lie group 
and Magnus expansion. In Section 3, we present the numerical procedure of 
Magnus expansion for solving matrix differential equations. In Section 4, we 
apply the Magnus expansion for computing eigenvalues of FSLP. In Section 
5, some numerical examples are given. 
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2 Magnus expansion 


The concept of Lie group was first proposed by Sophus Lie in 1893 to study 
transformation groups [13]. A Lie group (G,.) is a differentiable manifold 
that also has the algebraic structure of a group. For a given Lie group G, the 
tangent space at identity is called the Lie algebra of G and denoted by g. A 
binary operation [.,.] : gx g > g is called the Lie bracket. Let x,y € g, then 
the Lie bracket [z, y] is defined as 


[v,y] = xy — yz. (6) 


The operator (6) is bilinear, anti symmetric and satisfy the Jacobi identity, 


[x,y], 2] + [y, 2], 2] + [lz, 2], 9] = 0, (7) 


for all x,y,z € g. The Lie algebra g is a vector space under Lie bracket (6). 
An important special case of Lie group is the Special Lie group (SL(F,n)). 
Let F'”"*” be the set of all n x n matrices with entries in F’. The set 


SL(F,n) = {y € F"*", det y = 1}. (8) 


The usual matrix multiplication has the structure of a Lie group and is 
called the Special Linear group. The corresponding Lie algebra sl(F,n) of 
this group is the set of all n x n matrices with zero trace. The Lie bracket 
in this case is [A,B] = AB — BA, for A,B € sl(F,n). For more details 
see [4, 11,16, 22]. 

The Lie group is the differentiable manifold, therefor we can study or- 
dinary differential equations on a Lie group. A system of linear ordinary 
differential equations on SL(F,n) have the form [11, 22] 


y =G(a)y, tr(G(z)) =0 <> dety = 1. (9) 


Iserles and Norset [11] studied the solution of the linear matrix differential 
equation by the Lie group. In [22] the author introduce the collocation and 
Runge kutta type methods for the Lie group and present algebraic formulae 
for coefficients of these methods. Ros et al. [3], studied the Magnus expansion 
and some of its application in linear system of differential equations. The 
Magnus expansion represents the solution of (9) in the form 


y(x) = e? yo (10) 


where yo is the initial data of the problem (9) and 
x 1 x s 
o(a) = ih G(s)ds + 5 / i G(s), G(w)]duds 
P i/ / / [G(s), [G(u), G(v)||dvduds 
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4 af a [ 160). Gu), ow)]aeduas +. (11) 


In practice, we have to consider only a finite number of terms in expansion. 
The following theorem denotes the order of this approximation: 


Theorem 1. (Adapted from Theorem 6 in [11]) Let ||G(h)|| = O(h*), then 
P 
—Vaihy| = OCOD), (12) 


where 


o1= G(s)ds, 
0 


02 = ; [ [ie (s), G(u)]duds, 
ih f fe G(v)||dvduds 
a5 sf [ mn i OW) ldududs: 


Moreover, if & is an approximation of 0 such that ||o —a|| < O(h%) then 


a 


03 = 


lle — |] < O(n!). 


In this paper we suppose that p = 3, and approximate a(h) as follows 


m= [ Gio ih hee w)]duds 
ny, ff G(v)||dvduds 
dee ([G(s),G(u)], G(v)|dududs. (13) 
ahhh 


We explain the error of Magnus expansion for FSLP by the following Re- 
marks. 


Remark 1. By Theorem 1, if || G(h) ||= O(1) then the numerical Mag- 
nus series scheme produces an (at least) local order O(h°) approximation 
of o. A Magnus series integrator with local order p will have global order 
p —1 [12], thus Theorem 1 guarantees (at least) error of order O(h*) (if 


|| G(R) |= O(2)). 


Remark 2. In FSLP the coefficient matrix G is also a function of eigenvalue 
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A and for large eigenvalues, the norm of matrix G is : ||G|| = O(A). Thus for 
large eigenvalues, the error is of order O(\h*). 


3 Numerical procedure 


The first step for computing the solution y(x) = e7“) yo by Magnus series is 
the truncation of the series (11). As mentioned before we use only four terms 
(13). Next, we apply a numerical integration method for approximating the 
integrals (13). Here we use the two point Gaussian integration method. Thus 
we have 

J’ G(w)da = ®(G(cih) + G(ceh)), 


Jo fr[G(k), G@)ldédk ~ —¥Bn2[G(cih), G(eoh)], 

1 Se So G(s), [G(u), G(v)]]dvduds (14) 
$f" fe fIG(s), G(u)], G(v)|dvduds 

~ ([G(crh), G(coh)], G(erh) — G(coh)], 


Substituting (14) in (13) we obtain 


Oe *(G(erh) MOCK = vBip [G(cth), G(esh)] 
3 
+ £[1G(erh), G(c2h)], G(erh) ~ G(cah) (15) 


where c; = 5 v5, C2 = 5 | va are the nods of Gaussian quadrature, see [11]. 


Finally by Matlab software we compute y(h) ~ e?“yq. For computing y() 
in interval [a,b] we divide the interval [a,b] into n subintervals with length 
h= non and use the following iterative process 


Y (a;) = e?®9Y(a_1), Y(a)=Yo, 2 =a+in. (16) 


A Eigenvalues of fourth order Sturm-Liouville problem 


In this section, we formulate the fourth order Sturm-Liouville problem as the 
Lie group matrix differential equation and find the eigenvalues. For solving 
the system (4) and finding the eigenvalues of FSLP, we consider the matrix 
differential equation 
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Y’'=G(a)Y, Y(a)=Ih. (17) 


Theorem 2 [14]. Jf ®(x) is a solution of the matrix differential equation 
X'=G(a«)X on an interval J and if r is any point of J, then 


Vr,a2 € J, det ®(x) = det e(r)exp| tr(G(s))ds}. (18) 
From (5) it is clear that tr(G(x)) = 0, and by Theorem 2 for any solution 
Y (x) of system (17) we have 


det Y (x) = det V(expi [ tr(G(s))ds] = 1. 


a 


Thus system (17) is a system of matrix differential equations on the Lie 
Group SL(R*,4). Therefore we can solve the system (17) by Magnus expan- 
sion. Let Y (a) be a solution of (17) for x > a, then the solution of the system 
(4) with initial condition U(a) is 


U(x) =Y(2)U(a). (19) 


Thus, U(b) = Y(b)U(a). Following U(b) in the boundary condition (4) we 
obtain 
AU (a) + BY (b)U(a) =0 = [A+ BY (b)]U(a) =0. 


The initial condition U(a) is nonzero if U(a) = 0, then we have the trivial 
solution U = 0, thus determinant of the coefficient matrix A+ BY(b) must 
be zero. This determinant is a function of A and the eigenvalues of the FSLP 
are roots of equation: 


F(A) := det(A + BY(b)) =0. (20) 


Equation (20) is the main equation of this paper. In fact, this equation is 
the characteristic equation of the fourth order Sturm-Liouville problem. For 
solving system (17) and computing Y(b) by Magnus expansion, we need the 
brackets in (15). For FSLP these brackets are as follows 


[G(s), G(u)] = 
0 0 r(u)—r(s) 0 
0 r(s)qi(u)—r(u)qi(s) 0 r(s)—r(u) 
0 0 r(u)qu(s)—r(s)q1 (uv) 0 
0) 0 0) 0 


+ {A(w(u) — w(s)) + (d2(s) — ga(u))} os ie Fe 
0 
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[[G(s), G(@u)], G(s) — G(u)| = 


0 0 0 0 
0 0 2(r(s)—r(u))(r(s)ai(u)—r(u)ar(s)) 0 
0 2(a1(s)—a1 (u))(r(u)an (8) —7(s) a1 (u)) 0 0 
0 0 0 0 
0100 
+ (r(u) ~ r(s))(ar(s) — ar(u)) | 000° 
0001 
0000 
0000 
+ 2(r(s) — r(u){A(w(s) — w(u)) + (aa(u) ~ a2(s))} | 12°? 
0010 


For computing Y(«;) and finally Y(b) for finding the characteristic equa- 
tion (20), we should compute e7(“) which o(2;) is a matrix function of un- 
known parameter A. In general, we can not obtain e7(*i) explicitly. But for 
any given real number 4, we can compute it by Matlab or Maple software. 
For solving this problem, we limit the parameter A into interval [o, A*] and 
find the eigenvalues of the FSLP in this interval. For this purpose, we divide 
interval [Ap, A*] into m subintervals Ap < Ay < ... < Am = * with length 
hi = Avra and consider the following algorithm: 


1. i=0, 


i) 


. Compute Y(b) for \ = \; with iterative process (8), 


ww 


. Compute D; = det(A + BY(b)), if D; = 0, then A; is an eigenvalue, if 
D,Dj+41 > 0 go to step 4. If D; Dj+1 < 0, we have eigenvalue in interval 
[A;, Ai+1]. By applying Bisection method for function 
F(A) = det(A + BY (b)), we compute the eigenvalues in [Aj, :41]. 


4.i=i+1, ifi<m go to step 2, if i > m go to step 5, 
5. End. 


Remark 3. We can not obtain the explicit form of characteristic function 
F(A), but we can compute the values of F(X) for any given real number )j. 
On the other hand by (3) the roots of F(A) is simple hence the convergence of 
Bisection method is guaranteed. Thus in step 3 we use the Bisection method. 
Also, we can apply other numerical methods for computing the roots of F'()). 


Example 1. Consider the fourth-order Sturm-Liouville problem 
y — (0.02x7y')’ + (0.00014 — .02)y(x) = Ay(x), « € [0,5], (21) 


with the boundary conditions 
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y(0) =0=y""(0), 


It is easy to show that the eigenvalues of this problem are the square of 
the eigenvalues of second order Sturm-Liouville problem 


y(5) =0 = y'(5). (22) 


—u" (x) + 0.012?u(r) = pu(z), 
u(0) = 0 = u(5). 


x € [0,5], (23) 


The eigenvalues of problem (23) can be computed with Matslise [18]. We solve 
the problem (21) and (22) using Magnus expansion with n = 500 (h = 0.01) 
and approximate the eigenvalues in interval [0,650]. In Table 1 we listed the 
first 8 eigenvalues and compared them with other methods. 

In this example, for lower eigenvalues ||G|| = O(1), thus by Theorem 1 
and Remark 1, we have error of order (at least) O(h*). For large eigenvalues 
\|G|| = O(A) such that we must have error of order (at least) O(A\h*). The 
absolute errors in Table 1 confirm this order of error. If we compare the 
results with reference solution [18], it is obvious that our results is better 
than results of [1,2,5]. The Homotopy method [1], for eigenvalue \5 has 
large error. Also for eigenvalues Ag and A7 the errors on [2] are larger. 


Table 1: Eigenvalues and absolute errors Ad; of Example 1 

Ax Matslise [18] | AA, Our method | AA, [1] Arg [2] AAx [5] 
0.2150508644 2.9E — 11 44BF—9 |] 70FL—-11 | 28h —- 13 
2.7548099347 1.8£ — 10 L1IB-—6 | 83£-11 | 30h —- 12 
13.2153515406 2.1£ —11 44EF—-—5 | 5.8h-—-11 | 42h —-11 
40.9508197592 4.7E — 10 26£-3|)61fh-11 | 41B-7 
99.0534780635 5.2E —9 3.0L —1 | 758-8 
204.355732268 2.4E — 7 11.3 12H —3 
377.430420689 7.6E —6 12.73 
642.590868170 3.4E — 5 215.5 


Example 2. We consider problem (21) in Example 1 with the different 
boundary conditions 


y(0) =0=y'(0),  y(5) =0=y'(5). 


For the Magnus expansion result with n = 500, see Table 2. For this 
problem, we don’t have the exact eigenvalues, but since the equation of this 
example is the same equation of the Example 1, we expect the same error 
order of Example 1. Also there is an agreement between the results of our 
method and the other methods in Table 2. 


Eexample 3. 
beam: 


Consider the following Clamped-Clamped Euler-Bernoulli 
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Table 2: Eigenvalues of Example 2 


144.28062714 
280.60096225 
496.38280690 


144.28062804 
280.58602048 


144.28062684 
280.60096700 


Ap, Our method Ax [5] At [2] Ax [19] Ax [20] 
0.86690250 0.86690250 0.86690250 0.86690250 0.86690250 
6.35768645 6.35768645 6.35768645 6.35768645 6.35768645 
23.99274685 23.99274698 | 23.99274685 23.99274685 23.99274685 
64.97866759 64.97863592 | 64.97866760 | 64.97866759 64.97866761 


144.28062693 
280.60096374 


((1+2)2u"(2))” =A(1+2)F u(a), 
u(0) =0=w'(0), 


u(1) =0=u'( 


BR 
Ww 


(24) 


Problem (24) has the same spectrum as the following uniform Euler- 
Bernoulli beam problem: 


O<s<L, L=2(/2-1), 
o(L) =0='(L). 


(25) 


We called that the problems (24) and (25) are isospectral. For isospectral 
beam equation see chapter 12 in [7] and papers [6,21]. The eigenvalues of 
uniform beam (25) are the roots of the equation 


cos(8)cosh(8) -1=0, A= (Fy, BAD. 


Z (26) 


We compute the eigenvalues of problem (24) by Magnus expansion with n = 
300 and compare them with exact solution of equation (26) in Table 3. This 
problem has large eigenvalues such that ||G|| = O(A). Thus the errors are of 
order O(Ah*) and the results in Table 3 confirm the error order of Theorem 
1 and Remarks 1 and 2. 


Table 3: Eigenvalues and absolute errors AA, of Example 3 


k | Ap Our method Ap Exact AXk 

1 | 1062.777339630 | 1062.777339606 | 2.4F —8 
2 | 8075.518441201 | 8075.518441201 | 4.5F — 10 
3 | 31035.57009892 | 31035.57010021 | 1.3 —6 
4 | 84807.08316091 | 84807.08315850 | 2.4F —6 
5 | 189248.7474113 | 189248.7474331 2.2b —5 
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5 Conclusion 


In this paper, a method based on the Magnus Lie group method developed, 
which can be used to compute the eigenvalues of the FSLP. Numerical ex- 
amples were given to confirm the efficiency and accuracy of the method. 
Although the error grows for large eigenvalues, however the examples show 
that comparing to the other methods, our method has good accuracy for large 
eigenvalues. Example 3 shows that this algorithm is applicable for computing 
the eigenvalues of nonuniform Euler-Bernoulli beam equations. Ledoux [17], 
presented a modified Magnus method and obtained the large eigenvalues of 
second order Sturm- Liouville problem with error of order O($). It can be 
an open problem that one may modify the method of this paper in order to 
improve the results for higher index of eigenvalues. 
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